#!/usr/bin/perl
# 2012/11/19 Compare the peak locations

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [first bed file] -b [second bed file] -o [output pdf file]\n";
}

my $FILE_tmp1 = 'tmp1';
my $FILE_tmp2 = 'tmp2';


my %opt;
getopts("a:b:o:", \%opt);
my $FILE1 = $opt{a};
my $FILE2 = $opt{b};
#chr1	149450	149900	149690	24880.104	59.71	0	0
#chr1	365000	365390	365180	10971.672	61.641	0	0
#chr1	445850	446210	446030	8530.467	45.645	0	0
#chr1	456350	456770	456560	14966.053	48.099	0	0
#chr1	702920	703340	703130	11036.8	52.434	0	0
#chr1	717950	718370	718160	17702.709	61.471	0	0

my ($name1, $dir1, $ext1) = &fileparse($FILE1, '\..*');
my ($name2, $dir2, $ext2) = &fileparse($FILE2, '\..*');

my $FILE_out = $opt{o};

my @First = &readPeaks($FILE1);
my @Second = &readPeaks($FILE2);

sub readPeaks{
	my ($file) = @_;
	my @array;
	my $fh = IO::File->new($file) or die "cannot open $file: $!";
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/^chr/){
			my (@data) = split /\t/;
			my @temp = @data[0,1,2];
			push @array, \@temp;
		}
	}
	$fh->close();
	return @array;
}


&MinCalc($FILE_tmp1, \@First, \@Second, 'a');
&MinCalc($FILE_tmp2, \@Second, \@First, 'b');


sub MinCalc{
	my ($FILE, $array1, $array2, $header) = @_;
	my $fh_out = IO::File->new($FILE, 'w') or die "cannot write $FILE: $!";
	$fh_out->print("$header\n");
	foreach my $k(@$array1){
		my ($chr_r, $start_r, $end_r) = @$k;
		my $MIN = 3000;
		foreach my $l(@$array2){
			my ($chr_q, $start_q, $end_q) = @$l;
			if($chr_r ne $chr_q){
				next;
			}
			if($start_r <= $start_q and $start_q <= $end_r){
				$MIN = 0;
				last;
			}
			if($start_r <= $end_q and $end_q <= $end_r){
				$MIN = 0;
				last;
			}
			if($start_q <= $start_r and $start_r <= $end_q){
				$MIN = 0;
				last;
			}
			if($start_q < $start_r){
				my $distance = $start_r - $end_q;
				if($distance < $MIN){
					$MIN = $distance;
				}
			}else{
				my $distance = $start_q - $end_r;
				if($distance < $MIN){
					$MIN = $distance;
				}
			}
		}
		$fh_out->print("$MIN\n");
	}
	$fh_out->close();
}



my $Rcommand = <<END;
R --vanilla --slave <<EOF
	data1 <- read.table("$FILE_tmp1", header=T)
	data2 <- read.table("$FILE_tmp2", header=T)
	attach(data1)
	attach(data2)

	br <- seq(0,3000, by=50)

	pdf(file="$FILE_out")
	par(mfrow=c(2,1))
	hist(a, breaks=br, xlim=c(0,3000), xlab="Distance from $name1 to $name2", main=paste("total:", length(a), "peaks"), col='snow3')
	hist(b, breaks=br, xlim=c(0,3000), xlab="Distance from $name2 to $name1", main=paste("total:", length(b), "peaks"), col='snow3')

	dev.off()
EOF
END

 `$Rcommand`;
unlink $FILE_tmp1;
unlink $FILE_tmp2;
